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Abstract 

Interactions between surfaces and proteins occur in many vital processes and are crucial in biotechnology: the ability to con¬ 
trol specific interactions is essential in fields like biomaterials, biomedical implants and biosensors. In the latter case, biosensor 
sensitivity hinges on ligand proteins adsorbing on bioactive surfaces with a favorable orientation, exposing reaction sites to target 
molecules. Protein adsorption, being a free-energy-driven process, is difficult to study experimentally. This paper develops and 
evaluates a computational model to study electrostatic interactions of proteins and charged nanosurfaces, via the Poisson-Boltzmann 
equation. We extended the implicit-solvent model used in the open-source code PyGBe to include surfaces of imposed charge or 
potential. This code solves the boundary integral formulation of the Poisson-Boltzmann equation, discretized with surface elements. 
PyGBe has at its core a treecode-accelerated Krylov iterative solver, resulting in 0{N log N) scaling, with further acceleration on 
hardware via multi-threaded execution on opus. It computes solvation and surface free energies, providing a framework for studying 
the effect of electrostatics on adsorption. We then derived an analytical solution for a spherical charged surface interacting with 
a spherical molecule, then completed a grid-convergence study to build evidence on the correctness of our approach. The study 
showed the error decaying with the average area of the boundary elements, i.e., the method is (9(1/V), which is consistent with our 
previous verification studies using PyGBe. We also studied grid-convergence using a real molecular geometry (protein G B1D4'), 
in this case using Richardson extrapolation (in the absence of an analytical solution) and confirmed the (9(1/V) scaling in this case. 
PyGBe is open-source under an MIT license and is hosted under version control at https://github.com/barbagroup/pygbe In addi¬ 
tion, we prepared “reproducibility packages” to supplement this paper, consisting of running and post-processing scripts in Python 
to allow replication of the grid-convergence studies, all the way to generating the final plots, with a single command. 

Keywords: biomolecular electrostatics, protein surface interaction, implicit solvent, Poisson-Boltzmann, boundary element 
method, treecode. Python, CUBA 


1. Introduction 

Proteins interacting with solid surfaces fundamentally appear 
in many biological processes. Adsorption serves a key function 
in natural activities, like blood coagulation, and in biotechnolo¬ 
gies like tissue engineering, biomedical implants and biosen¬ 
sors. A full understanding of protein-surface interactions has 
remained elusive Glia, but adsorption mechanisms are gov¬ 
erned by surface energy and often the dominant effect is elec¬ 
trostatics. As a free-energy-driven process, protein-surface in¬ 
teraction is difficult to study experimentally |0, and thus sim¬ 
ulations offer a good alternative. Full atomistic molecular dy¬ 
namics simulations demand large amounts of computing effort, 
however, so we often must resort to other methods. 

Protein electrostatics can be studied via modeling approaches 
using the Poisson-Boltzmann equation and implicit-solvent rep¬ 
resentations. These models are popular for computing solva¬ 
tion energies in protein systems naia, but few studies have 
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included the effect of surfaces. Lenhoff and co-workers stud¬ 
ied surface-protein interactions using continuum models dis¬ 
cretized with boundary-element 012118] and finite-difference 
methods 0[IO], in the context of ion-exchange chromatogra¬ 
phy. They realized that van der Waals effects can be neglected 
for realistic molecular geometries and that the model is 
adequate as long as conformational changes in the protein are 
slight Em. 

The aim of this work is to develop and assess a computational 
model to simulate proteins near engineered surfaces of fixed 
charge, using implicit-solvent electrostatics. We have added 
the capability of modeling a protein near a charged surface to 
our code PyGBe, an open-source cod^that solves the Poisson- 
Boltzmann equations via an integral formulation, using a fast 
multipole algorithm and gpu hardware acceleration. Previously, 
we verified and validated PyGBe in its use to obtain solvation 
and binding energies, by comparing with analytical solutions of 
the equations and with results obtained using the well-known 
APES software In the present work, we derived an 
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analytical solution for a spherical molecule interacting with a 
spherical surface of prescribed charge, and used it to verify the 
code in its new application and study numerical convergence. 
Using the newly extended code, we also studied the interaction 
between protein G B 1 D 4 ' and a solid surface of imposed charge, 
and conduct a grid-convergence study using this more realistic 
surface geometry. 

We intend our new modeling tool to be useful in studying the 
behavior of proteins as they adsorb on surfaces that have been 
functionalized with self-assembled monolayers (sam), which 
are modeled within an implicit-solvent framework as surfaces 
with prescribed charge. One application is biosensing, where 
the target molecules are captured on the sensor via ligand 
molecules (for which antibodies are a common choice). Favor¬ 
able orientations of ligand molecules lead to greatly enhanced 
sensitivity of biosensors cma, because binding sites need to 
be physically accessible to the targets. Studies of protein ori¬ 
entation near charged surfaces might look at how orientation 
can be influenced by engineering decisions regarding surface 
preparation, to aid the design of better biosensors. We explore 
this application in a companion publication that obtains proba¬ 
bility of orientations for an antibody near a surface, as function 
of changing conditions on charge and ionic strength US). In 
this paper, we present the details of a new analytical solution 
for spherical charged surfaces and molecules, grid-convergence 
studies for the interaction free energy in this case, and grid- 
convergence studies for protein G B 1 D 4 ' alone and interacting 
with a charged surface. The detailed analysis of the model is 
complemented with a diligent effort for reproducibility and we 
deposit both input and results data in accessible and permanent 
archival storage, in addition to the open-source code. 



Figure 1 ; Sketch of the process for generating a 
solvent-excluded surface (ses): a protein molecule contains a 
set of atoms that define a radius upon applying a force field 
and a probe the size of a water molecule is rolled to define the 
SES. Qi is the protein region and 0.2 the solvent region. 




Figure 2 : Sketch of a molecule interacting with a surface; Qi 
is the protein, O2 the solvent region, F 1 is the ses and F2 a 
nanosurface with imposed charge or potential. 


2. Implicit-solvent model for proteins near charged sur¬ 
faces 


tential or charge. This new setup is sketched in Figure]^ and is 
described mathematically by the following equations: 


The implicit-solvent model uses continuum electrostatics to 
describe the mean-field potential in a molecular system. A typ¬ 
ical system consists of a protein in a solvent, defining two re¬ 
gions: inside and outside the protein, with an interface marked 
by the solvent-excluded surface (ses). The ses, beyond which a 
water molecule cannot penetrate into the protein, can be gener¬ 
ated by rolling a (virtual) spherical probe of the size of a water 
molecule around the protein (see Figure [^. Inside the protein, 
the domain has low permittivity (e = 2 to 4 ) and there are point 
charges located at the positions of the atoms. The solvent re¬ 
gion, representing water with salt, has a permittivity of e a; 80 . 
A system of partial differential equations models this situation, 
with a Poisson equation governing inside the protein and a lin¬ 
earized Poisson-Boltzmann equation governing in the solvent 
region. Appropriate interface conditions on the ses express the 
continuity of the potential and electric displacement, complet¬ 
ing the mathematical formulation. 

This model has been widely applied to investigate interac¬ 
tions between molecules, such as in protein-ligand binding. We 
are interested here in an extension of the model to consider in¬ 
teractions between proteins and surfaces with an imposed po¬ 


V^0i(r) = - ^ —5(r, r^) in solute (fli). 


in solvent (Q2), 
on interface F1, 


vV2(r) = K^4>2{r) 

01 = 02 
(901 502 

502 

02 = 00 or - 62-2— - ctq on surface F2, 

5n 


( 1 ) 


Here, 0 ,- is the potential corresponding to the region with 
permittivity e,, and 0o and ctq are the set potential or charge on 
the nanosurface. 


Boundary integral formulation. We express the system of 
partial-differential equations in o by the corresponding inte¬ 
gral equations along the interface and the nanosurface, F1 and 
F2. Many authors have used the boundary-integral representa¬ 
tion of the implicit-solvent model to compute solvation energies 
of proteins ifTTHTsHT^ 1^1^ 1 ^ 1131 . but apart from work led 
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by Lenhoff 0, we know of no studies that account for interact¬ 
ing surfaces in the system. 

Consider the setting in Figure]^ with prescribed potential at 
r2. The application of Green’s second identity on the first two 
equations of Q yields: 


01 +^“‘( 0 i,r,)- 

1 qk 

^1 ^ 47 r|rn^ - r^l 


on fii. 


02 - ^fy^( 02 ,r,) + j 

(1^02x2) = 0 onQ 2 , ( 2 ) 

where (piXj - 0,(rr,) is the potential in region Q, evaluated at 
the surface F^. K and V are dehned as 


^f/V(0f,r2) “ ^ ^ [GL/F(rn,, rr;)] 0;,r, dF, 

^ ^0/,r2GL/F(rnprr,)dF, (3) 

corresponding to the double- and single-layer potentials of 
and ^ 4 >iXj evaluated in the region Q,. The functions Gl 
and Gy are the free-space Green’s functions of the Poisson 
(Laplace kernel) and linearized Poisson-Boltzmann (Yukawa 
kernel) equations, respectively: 


GL(rni,rri) 


GF(rn 2 ,rr 2 ) 


1 

47r|rn, -rpj’ 
exp(-A-|rni - rrj) 
47 r|rn, - rp, | 


( 4 ) 


We then take the limits rn, —» rp^, —> rpj, —> rp^ on 

Equation (|^, and apply the boundary conditions: 4 '\X\ - 02 ,r,, 
= e2|i02,r, and ([>2x2 = 0o to get the following system 
of boundary equations: 


^ H- ii:[‘( 0 pp,)-y[‘ (^ 0 i,ri 
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- _ V 
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0i,ri 


tYi I ^ 


-^‘‘( 0 i.r.)+^nM^^ 




-/fy (0o) + Vy 1 — 02 X 21 - 0 onFi, 


fi ,,r2 / 9 


-Ky^i^iXi) + I 


ei 


^-/:^^( 0 o) + y;^(— 02 X 21 = 0 onF 2 . ( 5 ) 


Rearranging terms, we write Equation 0 in matrix form, as 
follows: 
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If the surface F2 has prescribed charge, corresponding to a Neu¬ 
mann boundary condition, -62^^02x2 “ '^’o, the equivalent 
derivation yields 
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^ii=0 47 r|rr, -rt| 
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( 7 ) 


The formulation detailed in this section differs from the work 
by Lenhoff and co-workers 013 because they consider an in- 
hnite charged surface, modeled using a modihed Green’s func¬ 
tion to account for the half-space domain. Lenhoff’s approach 
has the advantage that the charged surface does not require a 
mesh, but cannot be applied to any situation where the surface 
has non-planar geometry. An inhnite surface may not be a good 
model if the size of a device is comparable to the protein’s, like 
it happens with nano-structures. In that case, one needs to be 
able to represent the detailed geometry of a device via a surface 
mesh, and satisfy boundary conditions there, as detailed above. 

The system in (|3 can be extended to account for Stern lay¬ 
ers and solvent-hlled cavities by adding more surfaces or inter¬ 
faces. In our code, we deal with multiple surfaces in the manner 
presented by Altman and co-workers ED, as described in our 
previous paper 113 . 


3 . Methods 

5 . 7 . Discretization 

To numerically solve the system in Q, we discretize the 
boundaries into flat triangular panels and assume that cp and ^ 
are constant within those panels. The discretized form of the 
integral operators is as follows: 

4:d.sc (-^(rr)) = X ^ [Gdr, ry)] dF,, 

j=i 

N 

4-4.0 (|^<^(rr)) = Z |,0(rr2) £ Gp(r,, rp^dF,, ( 8 ) 
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where Np is the number of discretization elements on F, and 
0(rr,) and Jj<^(rr,) are the constant values of (p and on panel 
Fj (we are somewhat abusing the nomenclature here by reusing 
the symbol F, which previously referred to the complete sur¬ 
face). By collocating r, on the center of each panel, we get a 
linear system of equations that look just like (|^ or 0, but the 
coefficient matrix is formed by sub-matrices of size Np x Np 
rather than integral operators. Each element of a sub-matrix is 
an integral over one panel F^-, with r, located at the center of the 
collocation panel F;, as follows: 


Vi,,7= r Gi(rr„rr,)dF7. (9) 

The terms on the right-hand side and the unknown vectors in 
the discretized form of Equation 0 are sub-vectors of size Np. 
In this case, each element is the evaluation on the collocation 
panel F,, written as 




5n 
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47r|rr, - r*| 


5ii 
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i=0 
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47r|r; - r^r 


( 10 ) 


where r, is located at the center of panel F,. 

In our numerical solution, integrals are calculated in three 
possible ways, depending on how close the panel is to the col¬ 
location point. When the collocation point is inside the ele¬ 
ment being integrated, we use a semi-analytical technique, with 
Gauss points placed along the edges of the element f23[ p. 49, 
ff.l E^II . If the integrated element is closer than 2L from the col¬ 
location point —where L - ^2 ■ Aj for Aj the area of panel j — 
we use a fine Gauss quadrature rule, with 19 or more points per 
element. Beyond a distance of 2L, elements have only 1, 3, 4 
or 7 Gauss points, depending on the case. 


3.2. Treecode-accelerated boundary element method 

Most modern implementations of the boundary element 
method (bem) use Krylov methods to solve the linear system, 
usually a general minimal residual method (gmres), which is 
agnostic to the structure of the matrix. In practice, Krylov 
solvers for bem require 0{n ■ Np) operations to obtain the un¬ 
known vector, where n is the number of iterations to get a de¬ 
sired residual, and is much smaller than Np. The O(N^) scaling 
is given by a matrix-vector product (with a dense matrix) done 
in every iteration; this is the most time-consuming part of the 
algorithm, and makes bem prohibitive for more than a few thou¬ 
sand discretization elements. 

But when we inspect the approximation of the integrals in 0 
with Gauss quadrature rules, we see that the matrix-vector prod¬ 
uct has the form of an A^-body problem, similar to gravitational 
potential calculations in planetary systems. In this case, the 


Gauss quadrature points act analogously to planets (sources of 
mass) and the collocation points are analogous to the locations 
where the gravitational potential is computed (targets points). 
There are several ways to accelerate this kind of computations, 
for example fast-multipole methods ESll . treecodes ll26ll . and 
fast-Fourier-transform methods Il27ll . In our numerical solution 
(developed as the open-source code PyGBe), we accelerate the 
A^-body calculation with a treecode ll26l l28l . making this part 
of the algorithm scale as 6>(A^log A^) rather than 0{N^). 

The treecode algorithm groups the sources and targets in a 
tree-structured set of boxes and approximates interactions be¬ 
tween far-away boxes using a series expansion—a Taylor se¬ 
ries, in our case. This allows for controllable accuracy that de¬ 
pends on the number of terms used in the expansion and the 
multipole-acceptance criterion that defines the threshold where 
the distance between source and target is far enough to approx¬ 
imate the interactions with expansions. Details of our imple¬ 
mentation of the treecode in PyGBe can be found in our previous 
work ifT^ . 

3.3. Energy calculation 

Figure shows a system with three types of free energy: 
Coulombic energy from the point charges, surface energy due 
to F 2 and solvation energy. The Coulombic energy arises sim¬ 
ply from the Coulomb interactions of all point charges. This 
section describes how we compute the other two components 
of free energy in the boundary-element framework. 


Solvation free energy. When a protein is in a solvated state, 
surrounded by water molecules that have become polarized, its 
free energy differs from its state in vacuo by an amount known 
as the solvation energy. Its free energy again differs in the pres¬ 
ence of other structures in the solvent, e.g., other proteins or 
charged surfaces. In this work we use the term solvation en¬ 
ergy to more broadly mean the change in free energy of the 
protein from its state in a vacuum, to its state in the solvent 
with any other components or structures. In single-molecule 
settings, this definition of solvation energy coincides with the 
energy required to solvate the molecule. 

To calculate the solvation energy, the total minus the 
Coulomb potential is applied inside the protein, i.e.. 


Fsolv “ /-> I P (0total 0Coulomb) (11) 

^ Jn 

W, 

~ 'y ] qkiftotal ~ 0Coulomb)(f^)9 (12) 

/t=0 


where p is the charge distribution, consisting of point charges 
(which transforms the integral into a sum). The total 
minus Coulomb potential includes the reaction potential— 
representing the response of the solvent by polarization and re¬ 
arrangement of free ions—and any effects from the immersed 
surface. We can also interpret it as the potential generated by 
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the boundary F of the molecular region Q. Taking the first ex¬ 
pression of Equation Q and subtracting out the Coulombic ef¬ 
fect yields 

^reac.r, = (13) 

Equation ( [TT| ) requires evaluating 0reac for each point-charge 
location r^. We obtain this by discretizing Equation ( [T3| ) and 
using the solution of the linear system in Equation (|^ or Equa¬ 
tion 0 as inputs. 

Surface free energy. Chan and co-workers Il29ll^ derived the 
free energy for a surface with a set charge or potential. They 
describe the free energy on a surface as 



Eigure 3: Sketch of system solved with Legendre polynomials 
expansions. 


F - - GcCr^dT for set charge, and 

2 Jr 

F --- f Gpfgdr for set potential, (14) 

2 Jr 

where fo and (Tq are the prescribed potential and surface charge, 
respectively. The potential is given by f(cr, R, x) = GdR, x)cr 
for the first expression and the surface charge by cr{<p,R,x) = 
Gp(R, x)(p for the second one. This is valid because we are using 
a linearized Poisson-Boltzmann model. 

Using constant values of f and per panel, the discretized 
version of Equation ( [T4| takes the form 

1 

^ " 9 and 

i=i 

1 

7=1 

where Aj is the area of panel j, and cr = To obtain the 
surface free energy, we can plug in the solution of the system in 
Equation 0 or 0 to Equation 


4. Analytical solution 

It is possible to derive a closed-form expression for the free 
energy of interaction between a spherical molecule with a cen¬ 
tered charge and a spherical surface with imposed potential or 
charge, like the one sketched in Eigure There are such an¬ 
alytical expressions for interacting charged surfaces ED, and 
interacting spherical molecules with multiple point charges in¬ 
side ES, but not for a situation where surfaces and molecules 
interact. Having such an analytical solution is of great utility in 
the development of a computational model for protein-surface 
interaction, because it will allow for proper code verihcation. 


4.1. Expansion in Legendre polynomials 

The system of partial differential equations from Equation Q 
models the electrostatic potential field in the setting of Eigure 
Eollowing Carnie and co-workers JH], the axial symmetry 
lets us formulate the solution of Equation 0 as an expansion 
in Legendre polynomials: 


Interaction free energy. When there are two or more bodies in 
the solvent, they will interact electrostatically. In order to com¬ 
pute the energy of interaction, we need to take the difference 
between the total energy of the interacting system and the total 
energy of each isolated component, where the total free energy 
is given by 


Etotal “ Ecouiomb Eguj-face + EsoJv. 


(16) 


01 


^ c„r'i'P„(cos 6>i) H- 
«=0 


q 

4ne\ r\ 


on Qi, 


02 = ^a„k„(xri)E„(cos0i) 

rt-0 

oo 

-H ^ h„k„(xr2)P„(cos O2) on O2, 
«=0 


(18) 


The interaction free energy is 


_ ^assembly \ ^ ^comp, 

total / i total 


being P„ the n*-degree Legendre polynomial and the modi- 
hed spherical Bessel function of the second kind. 

(17) We make use of the following addition formula ll^ . 


where Nc is the number of components in the system and E™ j'’' 
is calculated over the isolated component i. 


k„(xr 2 )P„(COS 02) = ^(2m -H l)Bnmim(Kri)Pm(COS 0i), (19) 

m=0 
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to reformulate the expression for (^2 in Equation ([TSll as 


00 

02 =^a„^„(Afri)P„(cos6»i) 
n=0 

00 00 

+ '^b„ ^(2m + 1)B 

nm^m (Kri)P„,{cos 6 »i) 

n=0 m-0 

CO 

02 = ^ bnk„(Kr2)Pn(COS 62) 
n =0 

00 00 

+ ^ On ^(2m + l)B nm^m (j<r 2 )Pm(cos 02 ). (20) 

«=0 m=0 

Here, is the modified spherical Bessel function of the first 
kind; is defined by 


Bum — ^ ^ntn^n+m- 2 v{kR), (21) 

y=0 

where R is the center-to-center distance; and is given by 
the following expression, with r;^ = r(v) (in this context only) 
representing the gamma function: 


A'' 


r„,v 4 .o. 5 rn,-y+o. 5 rv+o. 5 (« + m- v)!(n + m-2v + 0.5) 
7 rr„,+„-y+i. 5 (n - v)!(m - v)!v! 


( 22 ) 

Legendre polynomials are orthogonal to each other, and 
is independent of 0. Thus, taking the inner product of the 
expressions in Equations ( [T 8 | ) and ( |20l l with Pj{cos 0,), where 
/ = 1 or 2 , yields 


Constant potential (p on T 2 . 

The application of the boundary condition on r 2 , 0 (r 2 ) = 0o, 
where 0 o is independent on 02 , gives 


OO 

^ a„( 2 ; + 1 )Bn ji j(Kd 2 ) + bnk„(Kd 2 ) 6 „ j = 0 odo j- (26) 

«=0 

Combining Equations ( |25| ) and ( |26] l yields the following system 
of equations for the coefficients a„ and b„ 


lA + LB = - 


ei q 
e .2 Ane\d^ 


e 


MA + IB = 0oe 


(27) 


where 


hn — ^jn 

An = a„ ^Kk'niKdi) - ^^^k„{Kdi) 
Bn bnknipd 2 } 


Ljn - (2y + 1)-B, 


Mjn = (2j+ l)B„jij(Kd 2 ) 


i'j(Kdi) ei j ij(Kdi) 
kn(Kd 2 ) e 2 d\kn{Kd 2 ) 
1 


(Kk'niKdi)- f^fkniKdl)) 


(28) 


Constant surface charge cr on Y 2 . 

In this case, the application of the boundary condition on r 2 , 
cr(r 2 ) = -e 2 g 5 lr 2 = fo, where ctq is independent on 6 * 2 , gives 


(pidQj = CjC. + — 6 oj (23) 

‘ 47reiri 

for the first expression of Equation ( [TS) , and 

OO 

(^2^0; ^ajkjiKrx) + ^ bn{ 2 j + 
n-O 

OO 

02^0; =bjkj{Kr 2 ) + ^ a„(2; + l)B„jij(Kr 2 ) (24) 


n=0 


for Equation ( |20l ). 

Applying the interface conditions for Ti on Equation 
and the first expression of Equation ( |24l i, produces 

'y anlKkniKdi) - — ^k„(Aft/i) 5„^ + 

\ > 

b„{2j + l)BnjiKij(Kdi) - --^ijiKdi) 


62 di 

e\ q 

62 An 6 id^ 


6 oj{j+l), (25) 


y an(2j + l)BnjKil{Kd2) + b„Kkn{Kd 2 ) 6 „j = -^<^07 (29) 


«=0 


Combining Equations ( |25] l and ( |26l l produces a system of equa¬ 
tions for the coefficients a„ and b„ 


IA-hLB = 

62 A7T6\d‘l 

MA-hIB = e 

e2 


where 


— ^jn 

An = a„ ^KkniKdl) - ^^^kniKdi) 
B„ = bnKk'n{Kd 2 ) 


Ljn - {2j + 1)B, 


nj 


i'jiKdi) ei j ij(t<di) 


yk'n{Kd2) 62 di Kk'n(Kd2) 
Mjn = (2j + \)BnjKlj{Kd 2 ) 


(30) 


[Kk'niKdy) - ||^k„(/^t/l)) 


( 31 ) 


where d\ is the radius of surface 1 . 
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4.2. Energy calculation 

Solvation free energy of the molecule. According to Equation 
o, the solvation free energy of a molecule with a centered 
charge is given by 

^’solv = ^g'0reac(ri = 0), (32) 

and using Equation ( [TS) , the reaction potential from Equation 
( pjj ) is: 

CO 

0reac = 0 “ = V P„(cos Ox). (33) 

Ane\ r ^ 

n=0 

Applying the boundary conditions at Ei on Equation ( |2^ , we 
can rewrite Cj in terms of the already computed aj and bf 

Cj = —(ajkj(Kdi)+ 

oo 

y, b„i2j + l)B„jij{Kdi) - ^ 6oj) (34) 

Because the charge is located at r = 0, only the n = 0 terms 
of Equation ( [33l l will survive, and the potential at this location 
is: 


<^reac(''l = 0) =ao^o('‘''^l) + 

OO 

y b„,B,„oio(Kdi) - ^ (35) 

The result from Equation ( |T5| ) in Equation ( |32| i yields the 
solvation free energy. 

Eor the isolated molecule, R oo makes B„m —» 0, which 
nullifies the sum in Equation ( |T5| ) and oq for —> oo, from the 
system in Equation ( |27l l, is 

PC ^ g gi 1 

d^ 62 4nKk'^(Kdi)ei 


Surface free energy with set potential (f>Q. We can expand Gp 
from Equation d in Legendre polynomials as 


Gp- 


62^ 

00 


y b„k'„(Kd 2 )P„(cos 02) 

n=0 


oo oo 

+ y a„ '^{2m + \)Bn„i„fKd2)Pmicas Oi) 

n=0 m=0 


(37) 


Applying Equation ([37]i in Equation (pA]) gives 


F - 2jiK<pod\e2 


bok'^iKd2) + y a„BnoiQiKd2) 


n=0 


(38) 


If the surface is isolated, R ^ oo makes B„o —» 0, and the free 
energy in this case is 


F = 2nK(l)odlb^k'oiKd2)e2 (39) 


where is taken from the system in (|27|l considering B„„j 
which results in 


, oo _ 00 

“ kQ(Kd2)' 


^ 0 , 

(40) 


Surface free energy with set charge ctq. We can expand Gc 
from Equation (pA]) in Legendre polynomials as 


Gc 


1 

O-Q 


y bnk„(Kd 2 )P„icos 02)+ 

n=Q 


y a„ y (2m + 1)B nm^m (Kd2)Pm(C0S 02) 

n=0 m=0 


(41) 


Applying Equation (|4T]l into Equation (p^ gives 


F - 2Kcr(jdl 


bokoiKd2) + y a„B„(iio(Kd2) 

n=0 


(42) 


Eor the isolated surface, R ^ oo and B„o —> 0, and the free 
energy is 

F = 2ncrodlb^koiKd2) (43) 

where is calculated from the system in ( |30l l considering 
Bnm —> 0, which results in 


b 


OO 

0 


tro 

e2KkQiKd2)' 


(44) 


5. Results of the grid-convergence study 

To obtain the following results, we extended the PyGBe code 
to consider surfaces with prescribed charge or potential. Eor all 
runs, we used a workstation with Intel Xeon X5650 cpus and 
one NVIDIA Tesla C2075 gpu card (2011 Eermi). We used the 
free msms software ll34l to generate meshes, and pdb2pqr llTSl 
with an amber force field to determine the charges and van der 
Waals radii. 

5.1. Verification against analytical solution 

Using the analytical solution detailed in Section]^ we car¬ 
ried out a grid-convergence study of PyGBe extended to treat 
interacting surfaces with biomolecules. The setup consists of a 
spherical molecule with a 5A radius and a centered charge of 
le , interacting with a spherical surface of 4A radius and an 
imposed potential of 0 = 1. The center-to-center distance be¬ 
tween the spheres is 12A, and they are dissolved in water with 
salt at 145mM, which gives a Debye length of 8 (ic = 0.125), 
and permittivity esoi = 80. The permittivity inside the spherical 
protein is fnioi = 4. Eigurej^ shows a sketch of this system. 

Figure [^presents the results of the grid-convergence analy¬ 
sis, where the error is the relative difference in interaction free 
energy between the analytical result from Section|^and the nu¬ 
merical solution computed with PyGBe. The observed order 
of convergence of the three finest meshes was 1.007. Table 
[T] presents the numerical parameters used in this case. Recall 
from sectionj^that we calculate the boundary-element integrals 
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Figure 1^ 


differently for close-by and far-away elements, and use a semi- 
analytical method for the element that contains the collocation 
point. The fine Gauss quadrature rule is used for elements 
closer than 2L from the collocation point, where L - V2 ■ Area. 
For the treecode, is the maximum number of boundary el¬ 
ements per box, P is the Taylor expansion truncation parameter 
and 6 is the multipole-acceptance criterion. The hnal numerical 
parameter is the exit tolerance of the gmres solver. 

Table 1: Numerical parameters used in the code-verihcation 
runs with the analytical solution. 


# Gauss points: 

Treecode: 


gmres: 

in-element close-by 

far-away A^crit P 

6 

tol. 

9 per side 37 

3 300 15 

0.5 

10-^ 



Number of elements 


Figure 6: Total runtime for obtaining interaction free energy 
(requiring three separate runs: one for each surface in isolation 
and one for both together), for the two-sphere system of|5.1| 
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Figure 7: Number of iterations to converge for the two-sphere 
system of 5.1 Iteration count increases with problem size due 
to the conditioning properties of the boundary integral 
formulation in our model. 



Number of elements 


Figure 5: Grid-convergence study for the interaction free 
energy between a spherical molecule with a centered charge 
and a sphere with potential (p — Data sets, figure hies plus 
running/plotting scripts are available under cc-by . 

As seen in Figure]^ the error decays with the average area of 
the boundary elements (2), which is the expected behavior and 
consistent with our previous verihcation exercises IfT^ . This 
proves that the extension of PyGBe to treat charged surfaces is 
solving the mathematical model correctly. 


Obtaining the interaction free energy involves three separate 
calculations: one with both bodies (molecule and interacting 
surface with set potential) and one for each isolated body. Fig¬ 
ure 1^ shows the total time to solution for each mesh, includ¬ 
ing all three cases that need to be computed. The most time- 
consuming part of the algorithm is the matrix-vector product 
within the Krylov solver, which scales as 0{N\ogN) thanks to 
the treecode acceleration. However, the total time-to-solution 
scales slightly worse than 0{N log N) because the condition 
number of the system increases with N, and we need more iter¬ 
ations of the Krylov solver to converge to the desired tolerance, 
as shown in Figure]^ 

5.2. Protein GBl D4' 

We computed the electrostatic field of protein GBl D4' in¬ 
teracting with a lOOAxlOOAxlOA block with surface charge 
density 0.05C/m^. The protein was centered with respect to a 
lOOAxlOOA face, a distance 2A above it. Since we did not 
consider any Stem layers or solvent-hlled cavities, these tests 
contain only two surfaces: the protein’s ses and the charged 
surface. We also computed the electrostatic field generated by 
protein GBl D4' and the surface by themselves. 
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^rot 


Table 2: Numerical parameters used in the convergence runs 
with protein GB1D4'. 


n2 



Figure 8: Orientation of a protein near a charged surface, m is 
the dipole moment vector, Vref the vector between m and the 
atom that is the furthest, and ni and n 2 are normal to their 
corresponding surfaces, atiit is the angle between ni and m, 
and (Trot the angle between Vref and n 2 when m and ni are 
aligned. 


The angle between the dipole moment of the protein and 
the vector normal to the surface was atilt = 10°. The dipole- 
moment vector placed at the center of mass of the protein gen¬ 
erates an axis, and we used the line of shortest distance between 
the outermost atom and this axis as a reference vector Vref- The 
rotation angle arot is the angle between the normal vector to a 
lOOAx lOA side face of the block and Vref when auit = 0, and is 
equal to 200° in these tests. Figure|^is a sketch of this arTange- 
ment. 

In this test, the solvent has no salt, i.e., k - 0, and its rel¬ 
ative permittivity was 80. The region inside the protein had 
a relative permittivity of 4. We computed the solvation and 
surface energy using meshes with 1, 2, 4 and 8 elements per 
square Angstrom with the parameters detailed in Table Us¬ 
ing Richardson extrapolation and the result of the three finest 
meshes we calculated an approximate exact solution, shown in 
Table|^ which we consider as a reference to calculate estimated 
errors. The errors plotted in Figures and are the relative 
differences between the energy obtained with Richardson ex¬ 
trapolation and the results computed with each mesh. They de¬ 
cay as 1/A in both the solvation and surface energies for the 
hnest three meshes, indicating that the calculations are in the 
asymptotic region and the geometry is well resolved in these 
cases. The observed order of convergence is 0.95 for the sol¬ 
vation energy and 1.12 for the surface energy in the isolated 
cases, and 0.96 for the solvation energy and 0.94 for the surface 
energy when the protein and surface were interacting. Using 
the extrapolated values from Table we obtain an interaction 
free energy of -7.6 [kcal/mol]. For details on the Richardson- 
extrapolation method for performing grid-convergence analy¬ 
sis, see our previous work El. 

5.3. Reproducibility and data management 

To facilitate the replication of our work, we consistently re¬ 
lease code and data associated with every publication. In that 
context, PyGBe was released under an MIT open-source license 


# Gauss points; 

in-element close-by 

Treecode; 

far-away A^crit ^ 

e 

GMRES; 

tol. 

9 per side 19 

1 500 15 

0.5 

10-** 


Table 3: Extrapolated values of energy for protein GBl D4'. 



Energy [kcal/mol] 


Solvation 

Surface 

Isolated 

-218.26 

317.41 

Interacting 

-222.43 

313.98 


with our previous publication IfTSlI . and continues to be avail¬ 
able via its version-control repository. Supplementing this pa¬ 
per, we prepared “reproducibility packages’’ containing run¬ 
ning and post-processing scripts in Python to generate Figures 
[^and[^ The packages invoke PyGBe with the parameters and 
meshes reported here, and then produce the plots, all with a 
single command. The reproducibility packages are hosted on 
figshare, and are referenced in the respective captions. 


6. Discussion 

In order to study the interaction of proteins and charged sur¬ 
faces, we extended PyGBe to account for surfaces with pre¬ 
scribed charge or potential. Unfortunately, there was no ana¬ 
lytical solution available in the literature to compare and verify 
PyGBe’s extension. Section [^derives a closed-form expression 
for a spherical molecule with a centered charge interacting with 
a spherical nanosurface with imposed charge or potential. We 
used this new analytical solution to conduct a grid-convergence 
study of the interaction energy (Figure |^. The error decays 
with the area, which is the expected behavior for a boundary el¬ 
ement method with constant elements Eiiia. Discretization 
error is very small for a spherical geometry. To make sure that 
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Figure 9: Grid-convergence study of the solvation energy for 
an isolated protein GBl D4' mutant, and the surface energy of 
an isolated surface with charge density of 0.05C/m^. 
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Figure 10: Grid-convergence study of the solvation and 
surface energy for protein G B1D4' mutant, interacting with a 
surface with a charge density of 0.05C/m^. Data sets, figure 
files and running/plotting scripts available under cc-by Eli 


Table 4: Numerical parameters for relaxed runs with protein 
GB1D4'. 


# Gauss points: 

in-element close-by 

Treecode; 

far-away A^crit ^ 

9 

GMRES: 

tol. 

9 per side 19 

1 300 4 

0.5 

10-= 


Table 5: Values of energy for protein GBl D4' using the 
parameters in Table and a mesh density of 4 elements per 
square angstrom in the protein and 2 elements per square 
angstrom on the charged surface 


Energy [kcal/mol] 
Solvation Surface 
Isolated -221.56 315.33 

Interacting -225.81 311.96 


the errors due to integration, the treecode approximation and 
the GMRES solver were even smaller, we chose all the numeri¬ 
cal parameters for high accuracy. This allows us to observe the 
convergence with respect to the discretization only and extract 
the order. With more realistic molecular geometries, however, 
discretization errors will be larger and the requirements can be 
relaxed in the other numerical parameters of PyGBe, resulting 
in lower runtimes. 


The results in Figures]^ and [T^ show the applicability of this 
approach in more realistic situations. The setting in Figure]^ 
can model a nanostructure coated with a self-assembled mono- 
layer (sam) interacting with a protein that will adsorb on that 
surface. Our application of interest in developing this model 
is the field of biosensors, where it can assist design through 
studies of electrostatic adsorption affecting protein orientation 
near the biosensor. With antibody-based biosensors, orienta¬ 
tion determines the accessibility of reaction sites and is critical 
for sensitivity. In a companion publication M, we present the 
first studies of protein orientation near charged surfaces using 
our modeling framework. 


Figures and 10 show the expected 1/N convergence with 
a simple protein. Just like in the case of the sphere, the sim¬ 
ple geometry of the charged surface forced us to use very fine 
parameters in order to extract the order of convergence. If we 
needed to run this computation many times—for example, to 
sample different protein orientations—we might relax these pa¬ 
rameters to obtain shorter computation times. 

The extrapolated values in Table are useful to find more 
relaxed parameters that still give acceptable results. For exam¬ 
ple, using a mesh density of 2 elements per square Angstrom on 
the charged surface and 4 elements per square Angstrom on the 
protein, and the parameters detailed in Table we get the re¬ 
sults in Table]^ These results are less than 2% away from those 
in Table and each run takes less than one minute. Moreover, 
using the energy values from Table|^ the interaction free energy 
is -7.61 [kcal/mol], which is very close to the extrapolated case 


7. Conclusion 

In this work, we used an implicit-solvent model to study 
protein-surface interaction. We present for the first time and 
apply an extension of our open-source PyGBe code to account 
for the presence of surfaces with imposed potential or charge. 
The new feature of the code was verified against an analytical 
solution, which we derived for that purpose. 

To demonstrate the power of this approach in a more realistic 
setting, we performed tests of protein GBl D4' near a brick¬ 
shaped surface with an imposed charge. The error in energy 
scaling with the area of boundary elements demonstrates that 
this extension of PyGBe is capable of resolving the mathemat¬ 
ical model correctly. This test was motivated by the biosens¬ 
ing application, where a ligand molecule is adsorbed on a sam- 
coated nanoparticle, which can be represented by the brick¬ 
shaped surface. 

The addition of a surface with imposed charge or potential in 
the implicit-solvent model falls naturally in a boundary integral 
approach. In this case, the region enclosed by the surface is not 
part of the domain, then, this surface only adds one equation to 
the linear system, rather than two, which is the case with the 
molecular solvent-excluded surface. 

We conclude that this implicit-solvent model can offer a valu¬ 
able approach in protein-surface interaction studies. This tool 
can be useful for orientation studies of ligand molecules in 
biosensors, either to find optimal adsorption conditions of salt 
concentration and surface charge, or to guide the design of bet¬ 
ter ligand molecules. 
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